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Using a phenomenological form of the equation of state of neutron matter near the saturation 
density which has been previously demonstrated to be a good characterization of quantum Monte 
Carlo simulations, we show that currently available neutron star mass and radius measurements 
provide a significant constraint on the equation of state of neutron matter. At higher densities we 
model the equation of state using polytropes and a quark matter model, and we show that our 
results do not change strongly upon variation of the lower boundary density where these polytropes 
begin. Neutron star observations offer an important constraint on a coefficient which is directly 
connected to the strength of the three-body force in neutron matter, and thus some theoretical 
models of the three-body may be ruled out by currently available astrophysical data. In addition, 
we obtain an estimate of the symmetry energy of nuclear matter and its slope that can be directly 
compared to the experiment and other theoretical calculations. 
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Introduction: While experimental information on mat- 
ter near the nuclear saturation density is plentiful, there 
are only a few experimental constraints on matter above 
the saturation density and, when available, they are con- 
taminated by strong systematic uncertainties. The vari- 
ation in the energy of nuclear matter with isospin asym- 
metry is particularly uncertain, since laboratory nuclei 
probe only nearly isospin-symmetric matter. There is a 
strong effort in trying to constrain the symmetry energy 
from intermediate energy heavy- ion collisions jlj, giant 
resonances in nuclei [2J, and parity- violating electron- 
nucleus scattering [3J 0] . 

Theoretical computations of neutron-rich matter are 
also difficult, owing to the poor quality of effective forces 
in dealing with neutron-rich matter [5] and uncertain- 
ties in the nature of the three-neutron force. At low 
densities, neutron matter is well understood because the 
two-body neutron-neutron interaction is constrained by 
experimental scattering phase shift data. At higher den- 
sities, three different classes of methods have emerged 
for computing the properties of neutron matter. The 
first class is based on phenomenological forces like the 
Skyrme interaction [B], The second class of calculations 
is based on microscopic nuclear Hamiltonians that typ- 
ically include two- and three-body forces obtained from 
chiral effective field theories, adjusted using renormaliza- 
tion group techniques to do perturbative calculations [7] . 
However the renormalization of the nuclear Hamiltonian 
induces many-body forces that have been carefully in- 
cluded in light nuclei [5] but not yet in nuclear matter. 
An alternative approach is the Brueckner-Hartree-Fock 
theory, which has been extensively used to study nu- 
clear matter and hyperonic matter [9lfTT| . The third class 
uses nuclear potentials, like Argonne and Urbana/Illinois 



forces, which reproduce two-body scattering and proper- 
ties of light nuclei with very high precision [TJJ [TS] . In 
the latter case, the interaction is designed to have small 
non-local terms, giving the potentials a hard core. The 
calculations can be performed in nonperturbative frame- 
work, and the strong correlations are solved by using cor- 
related wave functions. The ground state of nuclear sys- 
tems is determined by using the cluster-expansion |14j 
or using quantum Monte Carlo (QMC) methods. QMC 
methods have proven to be a very powerful tool to accu- 
rately study properties of light nuclei [TSl HH| and nuclear 
matter [T7] . All three of these classes suffer from strong 
uncertainties above the saturation density, both regard- 
ing the method and the nuclear Hamiltonian. 

On the other hand, astrophysical observations of neu- 
tron star masses and radii probe the equation of state 
(EOS) of dense, neutron-rich matter above the satura- 
tion density. Two types of neutron star mass and radius 
measurements have provided for progress on constrain- 
ing the EOS: the measurement of the general relativity- 
corrected radiation radius of quiescent low-mass x-ray 
binaries (qLMXBs) [TS], and the observation of photo- 
spheric radius expansion bursts which provides a simul- 
taneous measurement of both the mass and radius [T9l - 
[2T] . Reference [2U] has demonstrated that these two sets 
of data provide significant constraints on the EOS, rul- 
ing out several currently available theoretical models of 
dense matter. 

In this work, we show that these astrophysical ob- 
servations are beginning to constrain the nature of the 
three-body force in neutron matter. We construct a phe- 
nomenological description of the EOS near the saturation 
density which faithfully reproduces QMC simulations of 
neutron matter and can represent a wide range of EOSs 
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at high density [22, 23 . Utilizing the currently available 
astrophysical observations, we show that two parameters, 
closely connected with the strength of the three-body 
force and related to the magnitude and density depen- 
dence of the symmetry energy, are constrained by the 
observational data. 

The model: For densities below about half the satu- 
ration density, neutron stars consist of a crust which is 
solid except for a thin shell at the surface. Since the 
uncertainty in the EOS of the crust leads to an error in 
the radius which is much smaller than the current obser- 
vational uncertainty, we ignore variations in the EOS of 
the crust (see, e.g., [23]). We use the outer crust from 
Ref. [25] and the inner crust from Ref. [26]. Near and 
above the saturation density, we use a parametrization 
of neutron matter with the form 



£nm = P 



(1) 



where p is the nucleon number density, m n is the nu- 
cleon mass, and a, a, b and /3 are free parameters. In 
Ref. jTTJ, [22] it has been shown that this general form 
accurately fits the EOS of pure neutron matter given by 
QMC calculations using realistic nuclear Hamiltonians 
including two- and three-body forces. The uncertainty 
of the fit is much smaller than that from the three-body 
force. In the neutron matter case, the two parameters a 
and a are mostly related to the nucleon-nucleon force, 
while the parameter b is mostly sensitive to the corre- 
sponding symmetry energy E sym . The parameter /3 is 
sensitive to the particular model of three-neutron force 
(see Ref. [121 )• The range of parameters which subsumes 
all reasonable QMC calculations of neutron matter is 
12.7 < a < 13.3 MeV, 0.48 < a < 0.52, 1 < b < 5 
MeV and 2.1< j3 < 2.5. 

Because the parametrization in Eq. [I] describes neu- 
tron matter, we must also make a small correction to the 
neutron matter EOS due to the presence of a small num- 
ber of protons. In order to estimate this correction, we 
examine several Skyrme models from Ref. [57] , all chosen 
to have reasonable saturation properties and symmetry 
energies sufficiently strong as to prevent pure neutron 
matter from appearing in the maximum mass neutron 
star. We compute the mean and root-mean-square de- 
viation of the ratio of the pressure of neutron matter 
to the pressure of neutron star matter as a function of 
energy density over all the Skyrme models in our set. 
In our fiducial model, we apply a randomly distributed 
correction to the pressure with the same mean and root- 
mean-square deviation as that obtained in the Skyrme 
models, to our neutron matter EOS. This correction is 
0.86±0.03 at saturation density. For comparison, we also 
repeat our analysis without including the uncertainty in 
the ratio of the pressures, simply applying the mean cor- 
rection estimated from the Skyrme models. We have also 



checked this correction is similar to that given by similar 
relativistic mean-field models. 

Above the saturation density, four-body forces, hyper- 
ons, Bose condensates, and quark degrees of freedom may 
contribute to the EOS and our parametrization will no 
longer be appropriate. We take p t — 0.40 fm~ 3 as a rea- 
sonable upper limit for our parametrization of neutron- 
rich matter. Our fiducial model describes matter at 
higher densities by using a sequence of two piecewise- 
continuous polytropes, P = e 1 + 1 /™ with polytropic index 
n. The three parameters which describe the high-density 
EOS are m and ri2, the indices of the two polytropes 
and ep, the transition energy density between the two 
polytropes. Similar parametrizations have been used to 
describe the EOS at high densities and can mimic the 
presence of phase transitions. 

Alternatively, we describe matter at high densities with 
a parametrization of quark matter, with a polytrope at 
moderate densities to represent the possible presence of a 
mixed phase. For the quarks we use the model proposed 
by Alford et al. [2S]: 
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where P is the pressure, \i is the quark chemical poten- 
tial, the coefficient 0.6 < < 1 describes corrections 
to the massless free Fermi gas contribution from strong 
interactions, the coefficient a 2 — m 2 s — 4A 2 subsumes 
corrections from quark masses and color superconduc- 
tivity, and B is the bag constant. A largest possible 
range for a 2 is between (150 MeV) 2 -4(200 MeV) 2 which 
corresponds to a bare strange quark with a large quark 
gap and (400 MeV) 2 , which corresponds to a zero gap 
and strange quarks which receive significant contribu- 
tions from chiral symmetry breaking. This gives four 
parameters for the high-density part: the index of the 
polytrope, a 2 , 04, and the transition energy density be- 
tween the polytrope and quark matter which fixes the 
bag constant B. 

To match our parametrization to the constraints from 
neutron star mass and radius measurements, we use the 
method outlined in Ref. [20]. To that original data set 
we add a recent measurement from the transiently ac- 
creting neutron star U24 29J . We also add the constraint 
that models must be able to support at least a 1.93 so- 
lar mass neutron star, consistent with the 1 — a lower 
limit from Ref. [30]. We use Bayesian analysis, taking a 
uniform prior distribution for EOS parameters and using 
marginal estimation to compute the posterior probability 
distribution for EOS parameters, the EOS, and the mass 
versus radius curve. 

Results and discussion: Our constraints on the param- 
eters b and j3 for the neutron matter EOS near the satu- 
ration density are given in Fig. [T] and the parametriza- 
tion of the high-density part of the EOS is presented as 
Supplemental material [31| . We find that the posterior 
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FIG. 1. (Color online) The probability distributions and 68% 
(dark red areas) and 95% (green areas) confidence ranges for 
the parameters b and /3, and the density derivative of the 
symmetry energy, L. All distributions have been rescaled so 
that their peak is unity and then vertically shifted by an ar- 
bitrary amount. We compare our predicted value of L with 
constraints from nuclear masses ("Masses") [32], heavy ion 
collisions ("HIC") [l], pygmy dipole resonances ("PDR") [2], 
isobaric analog states in nuclei ("IAS") |33| . and antiprotonic 
atoms ("Pb(p,p)") |34| . The parametrization of the high- 
density part of the EOS is presented as Supplemental mate- 
rial En. 



probability distributions for the parameters a and a asso- 
ciated with the two-body force are almost flat as expected 
because they are related to the low-density part of the 
EOS. However, the parameters b and (3 are strongly con- 
strained by observations to ranges 3.3 < b < 4.8 MeV and 
2.28 < (3 < 2.5, independent of the nature of the high- 
density EOS. The results labeled "no corr. unc." contain 
no correction for the uncertainty in the ratio of the pres- 
sures of neutron matter to neutron star matter in each 
panel and are nearly indistinguishable from the results 
where this uncertainty is included (and as a result we al- 
ways include this uncertainty in the results shown below) . 
We have also checked that the effect of the varying the 
matching density the neutron matter EOS of Eq. [TJwith 
the polytrope between p t — 0.32 fm~ 3 and 0.48 fm~ 3 is 
small. 

Taking the nuclear matter binding energy at saturation 
to be —16 MeV, the nuclear symmetry energy is E sym = 
16 MeV + a + b. The density derivative of the symmetry 
energy, L = 3p (dE sym /dp) Po is 

L = 3(aa + b/3) (3) 

Taking the parameter ranges for a and b showed in Fig. 



[T]we find 32 < E sym < 34 MeV and 43 < L < 52 MeV 
to within 68% confidence. These tight constraints are 
possible because of an interplay between the strong cor- 
relation in QMC calculations of neutron matter between 
-E S ym and L and the constraints on the EOS from neu- 
tron star radii. The correlation between E sym and L re- 
sults from a separation between short- and long-distance 
parts of the three-neutron force [23] . Neutron star mass 
and radius measurements imply the radius is nearly in- 
dependent of mass and relatively small, which tends to 
select a smaller value for L and the value for E sym is 
then constrained from the correlation. We also find that 
the two- and three-body forces contribute almost equally 
to the symmetry energy at the saturation density, but 
this demarcation is more model dependent. Finally, our 
constraints on the neutron matter EOS are consistent 
with and complimentary to those from heavy-ion colli- 
sions [35) . which principally probe symmetric matter. 

The various models of three-body forces in neutron 
matter, that are typically constrained in light nuclei [13j . 
and the fact that the ranges of b and f3 are nearly indepen- 
dent of the high density EOS implies that neutron star 
mass and radius measurements can also constrain three- 
neutron forces. The ranges for b and j3 are smaller than 
the constraints determined from the wide range of pos- 
sible three-body forces given in Ref. [23], demonstrating 
that the astrophysical observations are ruling out more 
extreme models for the three-neutron force. The corre- 
sponding EOSs are given in Fig. [2] along with the EOS of 
Akmal, Pandharipande, and Ravenhall (APR) [14] and 
Skyrme model SLy4 [36] , The solid lines show the lim- 
its obtained by Gandolfi, Carlson, and Reddy (GCR) in 
Ref. [23] , obtained without any constraints from neutron 
star observations. 

In Fig. [3j we show the probability of the mass of neu- 
tron stars as a function of the radius for the different 
models. The quark models have slightly larger radii and 
slightly smaller maximum masses, but the general trend 
is similar to that described in Ref. [20 . Neutron star radii 
lie between about 11 and 12.3 km regardless of mass to 
within 68% confidence. As well as being consistent with 
QMC calculations from Ref. [53] these results are also 
consistent with the recent analysis of the EOS of neu- 
tron matter using chiral effective theories in Ref. [37 . 

Finally, we show the pressure of neutron star matter 
as a function of the energy density in Fig. [4] The results 
from APR [14], Sly4 [36] and GCR [23]. The quark mat- 
ter EOS is slightly softer at higher energy densities, but 
these energy densities are often beyond the central den- 
sity of the maximum mass neutron star. The pressure in 
APR is a bit larger than the observations suggest |20j . 
as is clear also in Fig. [2j We stress that any EOS out- 
side the results presented in Figs. [2] and |4j like APR at 
large energy densities, do not support constraints given 
by neutron star observations (see also Ref. [IS]). 

Conclusions: We find that neutron star mass and ra- 
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FIG. 2. (Color online) The 68% (dashed and dot-dashed 
lines) and 95% (shadowed areas) confidence ranges for the 
energy per baryon as a function of the baryon density as con- 
strained by the astrophysical observations. The APR [T^] and 
SLy4 [36] EOSs are also plotted, as well as the limits obtained 
in Ref. [23 (GCR). 
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FIG. 3. (Color online) The mass-radius curves for the mod- 
els considered in this work. The range of radii for 1.4 solar 
mass neutron stars, between 11 and 12 km, is similar to that 
obtained in Ref. [20] . The mass-radius curves for APR [14] 
and SLy4 [36] are also given. The labeling is the same as in 
Fig. @ 



FIG. 4. The pressure as a function of the energy density. The 
labeling is the same as in Fig. [2] 



dius measurements can be used to test and calibrate the 
three- and many-body nuclear forces in the contest of 
dense infinite matter. In particular we show that, when 
two- and three-body forces are parametrized as a sum 
of power-laws in the baryon density, the astrophysical 
observations strongly constrain both the coefficient and 
exponent which describes the three-body force. We also 
find novel constraints on the symmetry energy, driven 
partially by the strong correlation between S and L ob- 
tained in QMC calculations of neutron matter. 

There are potential corrections which are not yet well 
understood, including four-body forces, relativistic cor- 
rections, and the possible presence of hyperons [9TU1|. 
Our model partially takes these into account through 
the high-density polytropes (used above pt) which are 
not constrained. However, if these corrections are strong 
below n b = 0.32 fm~ 3 , then our constraints will have 
to be revisited accordingly. Some hyperonic models are 
consistent with our results [35] I3S] while others do not 
support neutron stars above 1.93 solar masses [101 [TT] . 

Another important difficulty is that there are several 
potential systematic uncertainties in the neutron star 
mass and radius measurements which are not yet under 
control. In the case of the photospheric radius expansion 
x-ray bursts, these include the nature of the relationship 
between the Eddington flux and the point at which the 
photosphere returns to the neutron star surface, the evo- 
lution of the spectrum during the tail of the burst [40] , a 
modification of the spectrum due to accretion, and viola- 
tions of spherical symmetry. In the qLMXBs, the X-ray 
spectra might contain high-energy power-law features not 
present in the atmosphere models. 
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SUPPLEMENTAL MATERIAL 



In this supplemental material section we provide the parametrization of the high-density part of the EOS described 
in the paper. The two different parametrizations, using two polytropes or one polytrope and a quark matter model, 
are presented in Figs. [5]and[6j 




FIG. 5. The probability distributions and 68% (dark red areas) and 95% (green areas) confidence ranges for parameters of 
the high-density EOS when represented by two polytropes. The polytropic index of the lower-density polytrope is m and 
the polytropic index of the higher-density polytrope is 712. The transition between these polytropes takes place at the energy 
density specified by ei. 




FIG. 6. The probability distributions and 68% (dark red areas) and 95% (green areas) confidence ranges for parameters of 
the high-density EOS when represented by a polytrope (with index ni) for the mixed phase and deconfined quark matter at 
higher densities. The transition between the polytrope and quark matter takes place at the energy density specified by e\. 
The parameter an and 02 represent the coefficients proportional to fi 4 and /i 2 in the pressure where fj, is the quark chemical 
potential. 



